////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//---------------------------- Ebrahim Foulaadvand- Somayyeh Belbasi, 13 Oct 2006 -----------------       //
//                                                                                                        //
//  The routine "Advection" solves the 1D wave equation using forward time centered space                 //
//  scheme and forward time backward space methods. Periodic boundary condition is used.
//  The system length L and the wave velocity c are set to one. Initial wave shape is
//  cosine-modulated Gaussian pulse.            //                                                         //
//                                                                                                        //
//                                                                                                        //                                                                                                        //
////////////////////////////////////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>

using namespace std;

main()
{

double tau=0.004,h,L=1,c=1,sigma=0.05,lambda=0.2,k=2*M_PI/lambda,x;

int N=200,i,n,T=1000;


ofstream file0 ("initial profile n=0.plt");  //output file for the profile at timestep n=0.
ofstream file1 ("profile FTCS n=20.plt");    //output file for the profile at timestep n=100.
ofstream file2 ("profile FTCS n=80.plt");    //output file for the profile at timestep n=500.
ofstream file3 ("profile FTBS n=20.plt");   //output file for the profile at timestep n=100.
ofstream file4 ("profile FTBS n=80.plt");   //output file for the profile at timestep n=500.
ofstream file5 ("profile FTBS n=100.plt");   //output file for the profile at timestep n=500.



vector<double> u(N+1,0),unew(N+1,0),uB(N+1,0),uBnew(N+1,0);

 // Arrays "u" and "unew" store the current and updated values
 // of the scaler function u(x,t) at grid points. "B" stands for backward space.

 h=L/N; cout<<"h= "<<h<<endl; //getch();

 cout<<"tau= "<<tau<<endl; //getch();


for(i=0;i<=N;i++){

x=-0.5*L+i*h; u[i]=cos(k*x)*exp(-0.5*(x/sigma)*(x/sigma)); uB[i]=cos(k*x)*exp(-0.5*(x/sigma)*(x/sigma));

file0<<(-0.5*L+i*h)<<"  "<<u[i]<<endl; }

//storing the initial condition data in the data file.



for (n=1;n<=T;n++){ //cout<<"n= "<<n<<endl; //getch();


// -------------------       FTCS method -------------------------------------

unew[0]=u[0]-(c*tau/(2*h))*( u[1]-u[N] ); // periodic condition has been used: u[-1] is replaced with u[N]

for(i=1;i<N;i++){ unew[i]=u[i]-(c*tau/(2*h))*( u[i+1]-u[i-1] ); }

unew[N]=u[N]-(c*tau/(2*h))*( u[0]-u[N-1] ); // periodic condition has been used: u[N+1] is replaced with u[0]


if (n==20)  for(i=0;i<=N;i++){ file1<<(-0.5*L+i*h)<<"  "<<unew[i]<<endl; } //storing in the data file.

if (n==80)  for(i=0;i<=N;i++){ file2<<(-0.5*L+i*h)<<"  "<<unew[i]<<endl; } //storing in the data file.


//-----------------------  FTBS Method      ------------------------------------

uBnew[0]=uB[0]-(c*tau/h)*( uB[0]-uB[N]);// periodic condition has been used: uB[-1] is replaced with uB[N]

for(i=1;i<=N;i++){ uBnew[i]=uB[i]-(c*tau/h)*( uB[i]-uB[i-1]); }


if (n==20)  for(i=0;i<=N;i++){ file3<<(-0.5*L+i*h)<<"  "<<uBnew[i]<<endl; } //storing in the data file.

if (n==80)  for(i=0;i<=N;i++){ file4<<(-0.5*L+i*h)<<"  "<<uBnew[i]<<endl; } //storing in the data file.

if (n==200)  for(i=0;i<=N;i++){ file5<<(-0.5*L+i*h)<<"  "<<uBnew[i]<<endl; } //storing in the data file.


for(i=0;i<=N;i++){ u[i]=unew[i]; uB[i]=uBnew[i]; }


}

cout<<"end." <<endl; getch();

}